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1. Accomplishments 


The status in the beginning of the report period was that the existing General Circulation 
Model (GCM) was running with a chemistry module compiled for stratospheric simulation 
studies. The chemistry simulation was not working sufficiently in the troposphere and any 
tropospheric trace gas sources or dry deposition sinks were not yet incorporated. 

The current status concerning the chemistry module is that 

• the chemistry simulation has been modified to also simulate the chemistry in the 
troposphere with resulting mixing ratios close to other model simulations as described 
in Olson et al. (1996). 

• The mechanism to incorporate trace gas source, e.g. testing for NO x , and dry 
deposition sinks, testing for H2O2, CH3OOH, O3, HCHO, HN 0 3) and NO2, are 
incorporated and is currently being tested. 

Existing model and development versions: 

• The full GCM model, currently still running with the original stratospheric chemistry 
module. 

• An off-line version of the GCM, i.e. wind and photolysis rates are pre-calculated and 
prescribed in read-in arrays. Here, the current modification are incorporated and are 
under test. 

• A box model version of the modified chemistry module for developments and first 
tests of new modifications. 

• A box model with same chemistry simulated but flexible partitioning and integration 
methods for test purposes of those. 

Several NASA scientist working in a joint effort on the development of the GCM at the 
NASA Langley Research Center which is funded through different grants. In particular 
three scientists, including the grantee, are working on the tropospheric chemistry part. 
Consequently the entire tropospheric development was not accomplished by the grantee 
alone. The grantee’s work focused mainly on three areas: 

• Trace gas sources and dry deposition 

• Method to introduce a NO source instead of a NO x source 

• Investigating integration methods 


These three areas and accomplishments are discussed in more detail: 



Tropospheric chemistry, dry deposition and trace gas sources 


A mechanism for computing trace gas production and loss terms caused by trace gas 
source and dry deposition fields has been incorporated in the off-line GCM version and is 
in a test phase. 

Mechanism: The calculation of deposition and source terms is using prescribed deposition 
velocity (currently for O3, HCHO, H2O2, CH 300 H, N 02 , and HNO3) and source flux 
arrays (currently for NO x ). Deposition velocities (or the associated transport resistance, 
see below) and source fluxes are the form in which , in general, the data are available in 
literature. 

The deposition velocity (v dep ) is defined as (Warneck 1990) 

Vdep = 1 / i^Rg “h Rs ) 

where Rg and R« are the ground and surface resistance, respectively, using the assumption 
that the surface flux is equal to the ground flux. R g reflects the transport resistance in the 
lowest few centimeters adjacent to a surface, and R« is surface specific, e g. surface could 
be bare ground or leaves. Deposition velocities are species and surface specific. In this test 
phase values from Muller (1992) for O3, HCHO, H2O2, CH 300 H, and HNO3 are used. 
The velocities for the different species are assigned to different vegetation classes. An 
longitude-latitude array has been generated using the NCAR vegetation database ds.769 
overlaid with snow and ice data used in the GCM. The resulting data array groups the 
earth surface into seven classes: ice/snow, water, bare ground, grass/shrubs, 

grass/shrubs/trees, non-tropical forest, and tropical forest (Figure 3) on a 0.5x0. 5 degree 
resolution. The deposition velocities are assigned to this database and then averaged to the 
128x64 model grid array. Resulting deposition arrays for HNO3, H2O2, and O3 are shown 
in Figure 4a-c. 

An additional transport resistance has to be introduced for the model surface interface, i.e. 
the column between surface and lowest model layer, in order to scale the deposition 
velocity which is defined at ground level (vaep) to a deposition velocity (v dep *) at the lowest 
model layer. This surface interface resistance Rj is the inverse of the model air/surface 
interaction coefficient D r (Model description; Stull 1988): 

Rj = l/Dr = U(Cd*M) 

where Cd and M are the drag coefficient and a mean wind velocity, respectively. V de p* is 
then (e g. Levy et al. 1985) 

Vdep* — 1 / ( Rg + Rd~) = Vdep X (1 — Vdep / Dr) 

D r , i.e. Cd and M, arrays will be provided from the transport routine of the model in order 
to calculated the modified deposition velocity Vdep* and associated loss rate in the 
chemistry routine. M is differently calculated for stable and unstable planetary boundary 
layers. Currently, for D r is assumed some global averaged number for the off-line test 
version of the model. 



The conversion of deposition velocities and sources fluxes into loss and production terms 
is performed as 

Ldep\ 1 / sec] = Vdep * !{dz / 2) 

Psrc{ molec / cm / sec] = F surface/ (dz/ 2 ) 

/Vc-f molec / cm / sec] = Fceii/dz 

where dz is the thickness of the model layer. For fluxes ad losses from/to the ground 
surface one has to use the column height between ground and the actual height of the layer 
center, that is dz/2, where the model calculations are performed. For sources generated 
inside a certain cell volume not originating from the surface, like the NO source from 
lightning or aircraft emissions, the thickness of the entire layer dz has to be used. As 
example, the difference between one-day model runs with and without deposition is shown 
for the resulting H2O2 and O3 mixing ratios in the lowest model layer in Figure 5. As 
expected from the deposition velocity arrays (compare Figure 4), the major losses appear 
over land and water and over land for H2O2 and O3, respectively. 

Status: The NO x source and deposition terms added to the off-line GCM version are still 
under test. 


• A method to introduce NO sources instead of NO x sources to the chemical routine. 

Problem: Since combustion processes emit mostly NO it would be more appropriate to 
introduce an NOx source in form of NO into a chemical model simulation. During 
daytime, NO is oxidized to N 02 and N 03 and rapidly recycled to NO through photolytic 
decomposition of NO2 and NO3 forming a catalytic chain reaction. Especially in the limit 
of nighttime conditions, i.e. under missing photolytic decomposition of NO3 to N02 and 
N02to NO, introducing NO x in form of NO2 instead of NO would miss the oxidation step 
NO + O3 and its associated loss rates for ozone. In the limit of daytime conditions it 
makes primarily no difference of adding either NO2 or NO because the catalytic NO-NO2 
cycle has a very high gain (chain length) in the order of several hundred recycles before 
NO/NO2 is lost into long lived reservoirs like HNO3, HN 0 4 , or PAN, and the difference of 
transitions N0-^N02 and N02~^N0 can only be maximal one of total hundreds per added 
NO or NO2. The transition between day and nighttime conditions is more complicated 
while the chain length of the NO/NO catalytic cycle drops to zero and the oxidation to 
NO3 and the loss to N2O5 increases. Nevertheless, since all this transitions are faster than 
tens of seconds and fast compared to integration time steps of typically dt> 900 s for global 
or regional model simulations, the NO-NO2-NO3 systems will be in photochemical 
equilibrium (PCE) and can be analyzed under fixed conditions, i.e. fixed mixing ratios 
resulting from the partitioning of NO x , for a particular time t. A NO source cannot simply 
be added as NO to the chemical system because NO-NO2 or NO2-NO3 already appear in 
PCE after partitioning and the fast relaxation processes, esp. in catalytic chains, cannot be 
resolved under a time step of dt. 



Method: The chemical system is analyzed for a certain time t, this is no time step t 0 ~^ti is 
performed but the flow of chemical reaction paths is analyzed based on the branching 
ratios between competing reactions under fixed species concentrations. The total sum of a 
species loss or production from reactions inside a feedback loops, like N0<^>N0 2 , can be 
parameterized in its loss or originating reaction times the gain (chain length) of the 
feedback loops. E.g. the loss of O 3 from its reaction with NO (LO 3 CNO+O 3 )) per initial 
NO molecule (O 3 yield per NO) can be written as the branching ratio of NO reacting with 
O 3 against the sum of all possible NO reactions (aNo+ 03 ) times the N0ON0 2 loop gain 
G 12 (indices 1 and 2 symbolize NO and N0 2 ; G 23 will be the loop gain of NO^NCh), see 
Equ. 2. The lower of the indices of LO 3 symbolizes its yield per initial molecule X, e.g. 
NO, and the upper one symbolizes the chemical system analyzed, e.g. here the NO<^ N0 2 

W 3 (NO + OJl m = G 12 x <x NO+ o, 

Equ. 2 O3 loss per initial NO from 
NO+O3 inside the N0<^N0 2 catalytic 
chain loop. 


L0 3 [N0 + 0-i)\ No ^ = X a jN0 2 X a NO+O t 

Equ. 1 O 3 loss per initial N0 2 from NO+O 3 
inside the N0<^>N0 2 catalytic chain loop. 

loop (no feedback from NO 3 included). Equ. 1 shows the O 3 loss of this reaction yielded 
per initial N0 2 . During daytime, (and for not too high O 3 concentrations) (Xno+o3 and ajK02 
are close to one (>0.99), thus LO 3 either per initial NO or N0 2 result in the same value, 
this is it makes primarily no difference adding NO or N0 2 during daytime. It turns out that 
the entire O 3 loss yield per initial molecule X (X=N0,N0 2 ,...) can be parameterized into 
the form of Equ. 3. F,g, and h are product terms of branching ratios a describing the path 
from the initial molecule to the target reaction (e.g. NO+O 3 or N 0 2 + 03 , see examples in 
Equ. 2and Equ. 1), and Gi 23 is an additional overall loop gain arising from the coupling of 
the two feedback loops N0ON0 2 and N0 2 <^N0 3 . 

LO i \ x = G ]2i x |G| 2 x f x \cc j j + G 23 x g x y 0 CjJ + G 12 x G 23 ^ h x \ct k y^ 


Equ. 3 Total O3 loss per initial X (X=N0,N0 2 ,..) from the N 0 < ^N 0 2 < ^N 03 

coupled catalytic loops. 



Model incorporation: The chemical routine (chemical solver) 
is initialized with mixing ratios of transported long lived 
families, i.e. groups of species, and the source flux arrays. 
The initial mixing ratio of a species family is then the result of 
chemical integration from the prior time step, plus the 
modification from transport calculations, plus an additional 
amount calculated from the source flux. The chemical solver 
can be seen separated in the partitioning of the long lived 
families into their shorter lived members, the estimation of 
chemical tendencies based on the chemical interaction of the 
partitioned species, i.e. computing of production and loss 
terms, and the final integration of the families (Figure 1). The 
partitioning step is in general iterated to gain stable solutions 
for, e.g. radicals like OH or NO while their production and 
loss terms depend strongly on each other. 

The new approach introduces an additional subroutine which 
analyses product yields of processes fast compared to time 
step dt, particular in this case, the 0 3 yield of the N0-N0 2 - 
NO3 catalytic cycle per NO or N0 2 added. This analysis 
would be placed after the partitioning (see Figure 1, shaded 
box). Here, the ratio RLO3 of the 0 3 loss yield per NO over 
the loss yield per NO x in PCE (i.e. NO/NO x =a, NO=aNO x , 
N02=(l-a)N0 x ) will be calculated (see Equ. 4b); NO x in PCE, 
this is the NO/NO2 PCE ratio for daytime, and N0 2 for 
nighttime. RLO3 can be used to correct the loss term LO3 of 
the chemical tendency for O3 (see Equ. 4c) which will be, as 
usual, in a first step estimated for dt based on total NO x in 
PCE (see Equ. 4a). I.e. the fraction of LO3 originating from 
the NO source will be multiplied by RLO3, gaining a final 
corrected chemical tendency for O3. 



Figure 1. 
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Equ. 4 

a) total N0 X = initial N0 X (from prior time step) + N0 X source. 

b) Ratio of O3 loss yield per NO over yield per NO x (in PCE) for the 
N0<=^N02 < ^N03 catalytic chain loops. 

c) Correction of the loss term of the 0 3 tendency which has been calculated for 
dt and NOx 104 which will correct the O3 loss yield for introducing a NO x 
source in form of NO. 

Status: To test Equ. 4.b, RLO3 has been calculated for mixing ratios of species compiled 
in box model runs initialized with 4 standard sets of species mixing ratios, altitude, and 
temperature: the 4 “IPCC cases” marine”, “land”, “plume-X”, and “free-troposphere” 
which were used for model inter-comparisons (Olson et al. 1996). Results are shown in 
Figure 6. As expected, RLO3 is equal to one for daytime. RLO3 significantly becomes >1 
when the NO3 photolysis rate drops below 10% of its noon value. At nighttime there is no 
photolytic decomposition of NO2 and NO3 recycling NO (■♦ G12— G23 = Gi23=l) Thus the 
maximum O3 loss per initial NO and NO2 is 2 and 1 (nominator and denominator of Equ. 
4.b), respectively, through reactions NO+O3 and NO2+O3. Since NO2+O3 is competing 
against N02+N03 - ^N 2 0 5 the nominator and denominator of Equ. 4b will range between 
1-2 and 0-1, resp. (neglecting NO+NO3), and RLO3 will become >2. The difference in 
amplitude of RLO3 between cases “land” and “marine” are due to different species mixing 
ratios, while the higher amplitude for “plume-X” and “free-troposphere” is due to lower 
temperatures and associated more favorable losses into the N 2 Oj reservoir. Increasing 
RLO3 values during night are associated with increasing NO3 mixing ratios (not shown) 
causing increasing losses towards N 2 0 5 . 

In a next step this needs to be added and tested in the off-line version of the GCM. 
Likely, this analytical method can also be used to calculate, e.g. OH, HO2, or O3 yields 
from fast non-methan hydrocarbons (NMHC) reaction chains when adding NMHC 
chemistry to the model. 



• Comparison of integration steps 


Integration in the chemical solver using 900 second time steps show variations after a 
compared to 10 second fine step integration, e g. after 5 days the differences for diurnal 
means can be in the order of 10% depending on initial conditions (IPCC cases). In the 
model version used for development, generally, for integration time step dt>5x (x is 
lifetime of species X), 0.2x<dt<5x, and dt<0.2x the integration over dt is calculated using 
photochemical equilibrium (PCE) assumption, implicit, and forward Euler integration 


X(t ] ) = (P-LxX(t 0 j)xdt 


*('.) = 


Pxdt + X(t 0 ) 


\ + Lxdt 

X(l,) = -*-(/„) + ( P / L - •*■('„)) x exp(- L X dl) 

Equ. 5 Integration steps: “forward Euler”, “implicit”, and “exponential”, respectively 


steps, respectively (see Equ. 5). All are functions of the production (P) and loss (L) terms. 
For a chemical solver as shown in Figure 1 with an iteration loop for the partitioning 
section the resulting terms P and L are generated for the time t 0 . The “correct” values for 
P and L would be some integral of P(t) and L(t) between to and ti which is unknown since 
calculations are only performed for times in step of dt. It has been investigated how 
integration results vary if one uses P(ti) and L(ti) instead of P(to) and L(to) or some 
average of them. P and L for ti can be calculated if one adds an additional iteration loop 
including the estimate of tendencies and the integration of families (compare Figure 1). 
The integration results for dt are compared to a fine step integration using forward Euler 
integration in steps of x/100 which don’t differ anymore from results using even finer time 
steps and are assumed to be “correct” solution with error «1%. Different P/L shapes 
were assumed between to and ti as illustrated in Figure 2. The envelopes of the resulting 
ratios between dt step integration and x/100 step integration are shown in Figure 7-9 as 
function of x. Positive numbers in these graphs mean, the dt integration overestimates 
values for X(ti), while negative numbers indicate an underestimation. It turns out that 
forward Euler and implicit integration steps for the above defined time domains using P 
and L of to general underestimate any results while using P and L of ti overestimate any 
results for X(ti) (compare Figure 7, Figure 8). This general trent could also be verified in 
box model runs using these two possible sets of P and L values. Using an exponential 
integration step (this is the correct solution for a P/L shape stepping at to to its final values 
P/L(ti) and staying constant over dt and only for this shape!) gives similar results as the 
implicit step with slightly smaller error for the implicit one (the results using the 
exponential step are not displayed). ( Note, that the large errors for the PCE values 
originate from concave shape calculations. A PCE state can hardly be reach if the change 



of the P/L value increases in time. At some point he system cannot anymore follow the 
changing P/L value fast enough. ) 
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Figure 2 Shapes of (production)/(loss)=P/L term steady state values 
assumed between integration time step to and tj. In “normal” the 
species mixing ratio at to is at its equilibrium value P/L(to), in “+/-“ 
cases the P/L value crosses the start value Xo. Shapes with opposite 
signs (mirrored at y=Xo will just flip sign of integration results. 
Concave and sinus-concave cases reflect the situation after sunrise 
and after noon where P/L values generally change increasingly. 
Convex and sinus-convex cases reflect cases befor noon and before 
sunset where changes in P/L decrease. Approximately linear changes 


Figure 9 show the envelope of results choosing P=P(ti) and L=L(ti) for the linear and 
convex P/L shapes, and P=(P(to)+P(ti))/ 2 and L=(L(to)+L(ti))/2 for the concave shapes, 
additionally implicit steps are used for dt>5x instead of PCE calculation in case of concave 
shapes. The error band indicated by the envelope curves is much smaller than for the 
calculation using either P and L only as function of to or tj. A more sophisticated 
exponential integration step, consisting of three different exponential terms which are 
function of P(t 0 ), P(ti),L(to), and L(ti) has been found which can compute results with 
much lesser errors than any other integration step discussed here. But two parameter have 
to be optimized for the specific shape of P/L. 


The point is, one can possibly improve the performance of the integration step by using 
not only production and loss term (P and L) as function of to but also of ti, as discussed 
above. But one has to introduce the mentioned outer iteration loop in the chemical solver 
routine which would mean that more computing time is used. Additionally, one has to 
setup criteria that define for what time period (associated with different P/L shapes in 
time) which kind of averaged value of P and L need to be used (two different types might 
be sufficient). This might be possible but has not yet been tested thoroughly with the box 
model, and improvements could not be verified yet. 









2. Comparison to proposed work 


The work proposed for the report period was mainly focused on preparing data arrays 
from field experiment, e g. O 3 LIDAR data, and starting a data-model comparison. The 
status of the tropospheric part of the model in the beginning of this report period was far 
from starting reasonable comparisons, e.g. the at this time incorporated basic chemistry 
was optimized for the stratosphere requirements and needed to be modified to work in the 
troposphere, neither were surface source and sink (dry deposition) terms incorporated 
which are essential for global tropospheric simulations. It seemed therefor to be more 
important first to focus on the model development to reach a status of the model where 
reasonable data-model comparison could be started. 

3. Milestones for the last grant year 08/96-07/97 

• Model development 

* Finalize testing of deposition and NO x source terms 

* Incorporate NO source instead of NO x source 

* Incorporate CO source and dry deposition, and also for hydrocarbons in time when 
associated chemistry will be incorporated. 

* Start to incorporate NMHC chemistry. 

* Look more into whether integration methods for chemical tendencies could be 
improved. 

* Develop and test above tools in off-line version and incorporate them in the full 
GCM version. 

• Data-model comparison 

* Prepare software tools and data arrays to use TOMS residual ozone database and 
the LaRC ozone sonde database for climatological model evaluations. Start model 
evaluations 

* Build up database of regional/large scale experiments, e g. TRACE-A/SAFARI- 
92 or PEM-West, for model comparison. 

* Use prescribed ECMWF or ECM wind fields with off-line GCM version for 
comparison with regional field data sets in specific time periods. 



Figure captions 


Figure 3. Vegetation classes grouped using the NCAR Olson ds.769 vegetation database, 
0.5x0 5 degrees (http://www.ucar.edu/dss/datasets/ds769-0. html) and snow/ice data used 
in the GCM. 

Figure 4, Deposition velocity grid arrays for a) HNO3, b) H2O2, and c) O3. 

Figure 5. Difference of a) H2O2 AND b) O3 mixing ratios in the lowest layer between 1 - 
day 3D model runs with and without deposition terms. 

Figure 6, Correction factor RLO3 for correcting the O3 loss term for adding a NO x source 
in form of NO instead of NO x . RLO3 is calculated following Equ. 4b for the four “IPCC 
cases”. Additionally, the NO2 and NO3 photolysis rates are displayed. 

Figure 7. Ratio of integration results as function of lifetime x of a species X. Used were 
P=P(tO) and L=(to). Integration performed in one time step dt=x over results in hundert 
time steps of t/100. Prescribed shape of P/L in time according to Figure 2 are used. The 
envelope of the results for every shape is shown, (V) for all shapes, (X) for shapes 
excluding “+/-“ cases and the concave cases. 

Figure 8. Same as Figure 7 but using P=P(ti) and L=L(ti). 

Figure 9. Same as Figure 7 but using P=(P(to)+P(ti))/2 and P=P(ti) for concave cases 
and for all others, respectively. (Same for L). For concave cases no PCE calculation is 
used but implicit steps for the entire range of dt>x/5. 
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